
# Read the execution times

mydata <- read.table("full.txt", header=TRUE, sep=" ")
nb.run <- dim(mydata)[1]
nb.fct <- dim(mydata)[2]

# Extract the empirical cumulative distribution function ECDF of the execution time of each user function

EmpiricalCDF <- NULL  
for (i in 1:nb.fct) {  
	EmpiricalCDF <- c(EmpiricalCDF, ecdf(mydata[,i]))  
}

ECDF <- function(i, x) {
	EmpiricalCDF[[i]](x)
}

# Plot all these ECDFs

for (i in 1:nb.fct) {
	plot(EmpiricalCDF[[i]],
		main = names(mydata)[i]
	)
}

#See Wikipedia
# Note that F_k^n(x) is EmpiricalCDF[[k]](x) because n is taken care of by ecdf

# C^n(u_1, u_2, ..., u_d) where u_k in [0,1]
Copula <- function(u) {
	total <- 0
	for (i in 1:nb.run) {
		result <- TRUE
		for (j in 1:nb.fct) {
			result <- result & (ECDF(j, mydata[i,j]) <= u[j])
			#print(paste("result = ", result, " and ECDF(", j, ",", mydata[i,j], ") = ", ECDF(j, mydata[i,j]), " and u[j] = ", u[j]))
		}
		if (result) total <- total + 1
		#print(paste("total = ", total))
	}
	total/nb.run
}

# Joint distribution F(x_1, x_2, ..., x_d) where x_k is the execution time of the k'th function
JointCDF <- function(v) {
	values <- NULL
	for (i in 1:nb.fct) {
		values <- c(values, ECDF(i, v[i]))
	}
	#print(paste("values = ", values))
	Copula(values)
}


#======================================================================#
#======================================================================#
#======================================================================#
#======================================================================#

#See the pdf www.eurandom.tue.nl/events/workshops/2011/MRMlectureday/Presentations/Can190111.pdf, slide 12/40
# Compute the inverse of the ECDFs
InvECDF <- function(j, p) {
	minimum <- -1
	for (i in 1:nb.run) {
		if (ECDF(j, mydata[i,j]) >= p) {
			if (mydata[i,j] < minimum || minimum == -1) minimum <- mydata[i,j]
		}
	}
	minimum
}

Copula2 <- function(p) {
# compute the inverse
	#print(paste("p = ", p))
	values <- NULL
	for (j in 1:nb.fct) {
		values <- c(values, InvECDF(j, p[j]))
	}
	#print(paste("values = ", values))
	total <- 0
	for (i in 1:nb.run) {
		result <- TRUE
		for (j in 1:nb.fct) {
			result <- result & (mydata[i,j] <= values[j])
		}
		if (result) total <- total + 1
	}
	total/nb.run
}

vecD <- c(21, 27, 145, 2110)
vecP <- c(0.1, 0.9, 0.3, 0.25)


#======================================================================#
#======================================================================#
#======================================================================#
#======================================================================#

#mydata <- read.table("output_10_100.txt", header=TRUE, sep=" ")
mydata <- read.table("full.txt", header=TRUE, sep=" ")

nb.run <- dim(mydata)[1]
nb.fct <- dim(mydata)[2]

observed.WCET <- rowSums(mydata)
observed.WCET.ecdf <- ecdf(observed.WCET)
plot(observed.WCET.ecdf)

sample.step <- 0.2
values <- seq(0,1,sample.step)
sample.point <- permn(values)
#sample.points <- matrix(runif(nb.point * nb.fct), nb.point, nb.fct)
estimated.copula.ecdf <- C.n(u=sample.points, U=pobs(mydata))
 
#EX <- 0
#x1 <- C.n(u=mat, U=pobs(mydata))
#x2 <- array(0, dim=c(m))
#for (i in 1:m) {
#	x2[i] <- Copula2(mat[i,]) 
#	if (x2[i] != x1[i]) {
#		print(paste("x1=", x1[i], ", x2=", x2[i]))
#		EX <- mat[i,]
#		print(EX)
#	}
#}


#y <- dCn(u=mat, U=pobs(mydata), j.ind=1:nb.fct, b=1/sqrt(nrow(pobs(mydata))))